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Abstract 

The calculation of quantum canonical time correlation functions is considered in this paper. 
Transport properties, such as diffusion and reaction rate coefficients, can be determined from time 
integrals of these correlation functions. Approximate, quantum-classical expressions for correlation 
functions, which are amenable to simulation, are derived. These expressions incorporate the full 
quantum equilibrium structure of the system but approximate the dynamics by quantum-classical 
evolution where a quantum subsystem is coupled to a classical environment. The main feature of the 
formulation is the use of a mapping basis where the subsystem quantum states are represented by 
fictitious harmonic oscillator states. This leads to a full phase space representation of the dynamics 
that can be simulated without appeal to surface-hopping methods. The results in this paper form 
the basis for new simulation algorithms for the computation of quantum transport properties of large 
many-body systems. 
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1 Introduction 

In the investigation of condensed phase systems one is usually interested in the average value of some 
observable or in a time correlation function from which a transport coefficient can be computed. Since 
the basic description of matter is quantum mechanical, we are interested in the average value of a quan- 
tum mechanical operator, which is given by B{t) = TiB{t)p{0), where, for a system with Hamiltonian 
H, the time dependent operator satishes the Heisenberg equation of motion, 

^m = liH,B], (1) 

and p{0) is the initial value of the density matrix. Quantum time correlation functions of two operators, 
A and B have the fornix, CAB{t) = Tr(/5cqAB(t)), where peq is the quantum canonical equilibrium 
density matrix. Either of these quantum expressions requires a knowledge of the time evolution of 
a quantum operator in a many-body system that is often very large. Consequently, these general 
expressions are not computationally tractable and appeal must be made to approximations if they are 
to be evaluated for problems of physical interest. 

The approximate dynamical description we consider in this paper is quantum-classical Liouville 
dynamics [2]. In this formulation the system is partitioned into two subsystems, which we call quantum 



^ Other forms of quantum correlation functions are useful in applications. These include symmetrized and Kubo 
transformed forms. Since relations exist among these correlation functions [1], we restrict our considerations to this 
expression. 
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subsystem and bath or environment. The partition is dictated by physical principles. For example, in 
electron or proton transfer problems, the electron or proton may constitute the quantum subsystem 
while the environment in which the transfer takes place, a molecular group or biomolecule dissolved 
in a solvent, forms the bath. In quantum-classical Liouville dynamics the bath is treated classically 
while retaining the full quantum character of the quantum subsystem. In this theory the analog of the 
Heisenberg equation of motion for an operator is [21 [3] 

—Bw{X,t) = - [Hw,Bw] - -{{Hw^Bw} - {Bw.Hw})-, (2) 

where X = (R, P) are the positions and momenta of the bath degrees of freedom. The subscript W 
indicates that operators are expressed in a partial Wigner representation defined below. In addition 
to the usual quantum commutator, the equation of motion also involves a Poissson bracket denoted 
by {•, •}. While this equation is far more tractable than the full quantum equations of motion, its 
numerical simulation presents challenges and a number of different schemes have been constructed for 
this purpose [H [5]. These include Trotter-based methods that use an adiabatic basis [6], as well as 
trajectory methods that use the diabatic and force bases [8], schemes based on the multiple threads 
algorithm [9] and a method that utilizes the mapping basis [10]. The mapping method has proved to be 
especially promising for the evaluation of expectation values of operators, as shown by calculations of 
population relaxation in the spin-boson model, one of the standard test cases for quantum dynamics [10] . 
In this paper, we show how the mapping method can be extended to the computation of time correlation 
functions within the quantum-classical Liouville framework. 

In the next section, we present the explicit expression for the quantum correlation function that 
forms the basis of our calculations. Section [3] outlines the mapping formulation where subsystem 
quantum states are replaced by fictitious harmonic oscillator states. Section ^reformulates the quantum 
correlation function in the mapping basis. This fully quantum description is exact and the dynamics is 
embodied in the spectral density function. A quantum-classical approximation for the spectral density 
dynamics is derived in Sec. El The details of this derivation are presented in the Appendix. This result 
is used in Sec. [6] to obtain the final result for the quantum-classical correlation function in the mapping 
basis. The conclusions of the study are presented in Sec. El 



2 Correlation Function 

Before considering quantum-classical approximations to the dynamics, we first rewrite the exact quan- 
tum correlation function in a form that is suitable for the introduction of the mapping basis and passage 
to the quantum-classical limit. Suppose the quantum subsystem and bath have Ng and Nf, degrees of 
freedom with characteristic masses m and M, respectively. In order to introduce a phase space de- 
scription of the bath, we first introduce a coordinate representation of the bath so that the correlation 
function takes the form, 

CABit) = Tr(peqAB(t)) = Tr(peqAe'"*/^Be-^^*/^) 

= Tr'y dQidQ2dQ3dQ4(Qi|/5eqAIQ2)(Q2|e*^*/'"^|Q3)(Q3|BlQ4)(Q4|e-'^*/f^|Qi) 

= Tr' I dRidR2dZidZ2(Ri - y|/OeqA|Ri + ^)(Ri + ^le'HV&iR^ _ ^) 

x(i?2 - ^\B\R2 + ^){R. + f |e-^*/^|i?i - f ). (3) 

In these equations, Tr' stands for a trace over the quantum subsystem degrees of freedom and, in the 
forth equality above, we have made a change of variables Qi = Ri — Zi/2, Q2 = Ri + etc. 
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Two further manipulations are required to cast the correlation function into a form that is suitable 
for the calculations using the mapping basis described below. The coordinate space matrix elements may 
be replaced with bath phase space functions by introducing the partial [3] Wigner transforms [IH [T2l [T3] 
of an operator and density matrix, 

(-R-flft,/i|fl + |) = j dP{p,,A)w{R.P)e--''-^l''. (4) 

Note that the partially Wigner transformed quantities are still operators in the quantum subsystem 
Hilbert space. In addition, Tr' may be written explicitly using the eigenfunctions of the Hamiltonian 
of the quantum subsystem. The Hamiltonian for the entire system is given by 

IM 2m 

where P and p are momentum operators of the bath and subsystem and Vf^, Vg and Vc are, respectively, 
the bath, subsystem and coupling potentials. The coordinate operators for the subsystem and bath are 
q and R, respectively. Equation ([5]) can be written as H = ^ + Vb + Vc + hs, where hs = + Vg is 
the subsystem Hamiltonian. The eigenstates of are defined by the eigenvalue problem hs\X) = e^lA). 
We suppose that there are A'' quantum subsystem states. Making use of subsystem energy eigenstates 
as a basis and introducing the Wigner transformed forms of Peq^ and 13 given in Eq. ^ we have, 

1 ^ r 

CAB{t) = (^^^ J dXidX2dZidZ2 {X\{pegA)w{Xi)\X') 



X,X',fJ.,fl': 

Z, . . Zo 



X 



The spectral density function in the subsystem basis may be defined as 
VF^'VM(Xi,X2,t) = j dZidZ2 {X'\{Ri + ^|e^^*/^|i?2 - 

xif^'m + f |e-^*/^|i?i - ^)\X)e-rAP-'^^P-^^\ (7) 

and contains all information needed to compute the quantum time evolution contribution to the cor- 
relation function. The spectral density defined via full Wigner transform was used by Filinov et al. 
|141 [TSl [IB] in a reformulation of the quantum correlation function. In terms of the spectral density 
function in the subsystem basis the correlation function takes the form, 

1 ^ r 

CABit) = j^—rj^ Yl / (A|a,i)H/(Xl)|A')(^|^H.(X2)|/x')H^^'^'^''^(^l,X2,t). (8) 

This expression for the quantum correlation function is exact. In order to compute it, the matrix 
elements of the forward and backward propagators must be evaluated to solve for the time dependence 
of the spectral density function. This is the most difficult part of the problem. A similar expression 
utilizing the adiabatic basis in place of the subsystem basis was derived earlier |17j . In addition to 
the time dependent spectral density function, the time independent matrix elements of the quantum 
operators and quantum equilibrium density matrix must also be computed to evaluate the ensemble 
average appearing in the definition of the correlation function. While these equilibrium quantities may 
be difficult to evaluate for complex systems, they are far easier to compute than the quantum time 
dependence, and algorithms have been developed for this purpose [19]. 



3 



3 Mapping Basis 



In order to construct a useful simulation algorithm for the quantum correlation function, it is conve- 
nient to use an alternative, but equivalent, mapping form for the quantum subsystem matrix elements 
of the operators which enter in its definition. A well-known mapping approach was introduced by 
Schwinger [20j. In his scheme, the eigenstates of the angular momentum operator are mapped onto 
eigenfunctions of two bosonic oscillators, and the angular momentum operators are mapped onto com- 
binations of creation, a\, and annihilation, oa', operators (A, A' = 1 or 2). Such a mapping yields a 
simple treatment of the angular momentum problem in quantum mechanics. In this formalism the 
resolution of identity is mapped to a\ai + 0202 = 1- This equality is used in the Holstein-Primakoff 
mapping scheme to eliminate one bosonic oscillator and represent angular momentum states by a single 
oscillator [21] . An extension of these mapping schemes was used to map discreet states of a quantum 
system onto fictitious harmonic oscillators, so that all degrees of freedom in the system can be treated 
with semiclassical approximations [22l [23l [24l ESI [26l [271 121] • The mapping approach has also been 
used in other closely related contexts [291 EHl [31] . 

In the mapping scheme used in this work, an N level quantum system is mapped onto N harmonic 
oscillators. The wavefunction of the system in a given quantum state is mapped onto a product 
of harmonic oscillator wavefunctions where all oscillators are in their ground state, except for one 
oscillator corresponding to the given quantum state, which is in its first excited state. This mapping 
is schematically represented in Fig [TJ Therefore, we have a physical space with a cardinality N, which 
is much lower than the cardinality of the Hilbert space of the harmonic oscillators, which is infinite. 
More specifically, the subsystem quantum states are mapped through the relations 

|A) ^ |mA) = |0i,--- ,1a,--- ,0^), (9) 

where 

{q\mx) = (91,92,- •• ,qn\Oi,--- ,1a,--- ,0n) = Mqi) ■ ■ ■ 4>o{qx^i)(l)i{qx) ■ ■ ■ (po{qN), 

with (pQ and (pi, respectively, being the ground and the first excited state wavefunctions of an harmonic 
oscillator. The creation and annihilation operators on the mapping states act in the following way^: 



aj^lO) = |0i,--- ,1a,--- ,OAr) = |mA), and OA|mA) = |0i, - - - , Oat) = |0), (10) 



where 



ax = \^{qx H Px), a{ = \Hz^{qx Px), and [qx,Px] = ifi- (H) 



We may then introduce the mapping representation of an operator A as 

= ^^AA'4"A'- (12) 

AA' 

For example, the mapping form of the Hamiltonian in Eq. ([5]) is 



P 

AA' 

72 



^ + V^B(i?) + ^ E ^AA'(^) (rxrx' + ^ - —5xx), (13) 

XX' 



■^Since the harmonic osciUators are fictitious we may set muj to 1 to simphfy these expressions; however, we retain these 
forms both to malte dimensional consistency manifest and to simplify the dynamical relations (see Sec. [SJ by equating m 
to the quantum subsystem characteristic mass and u) to the inverse of the scaling time. 
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0,0,0,1,0) 



Figure 1: Schematic representation of the mapping for a 5 level system. Subsystem states are on the 
left and the mapping states are on the right. 

where h = + Vs{q) + Vc{R,q) and /iaa' is its matrix element. We used the fact that hw = hyx 
in writing this expression. Note that unlike Eq. ([5]) where q and p are the subsystem coordinates and 
momenta, in Eq. (|13p r and p are the mapping space coordinates and momenta. These quantities have 
dimensions equal to the number of subsystem quantum states. 

From the definition in Eq. (jl2p it is clear that the matrix elements of A in the subsystem basis 
are identical to those of Am in the mapping basis: (A|A|A') = Axy = {mx\Am\mxi) ■ Consequently, 
any matrix element in the quantum subsystem basis can be substituted by its equivalent form in the 
mapping basis. We show that this substitution leads to computational advantages when simulating 
quantum correlation functions. 

4 Correlation Function in Mapping Basis 

The matrix elements in the correlation function expression ([6j) can be replaced by their mapping 
equivalent forms to yield 



The next step in the analysis of this correlation function is to introduce a coordinate space representation 
of the abstract mapping eigenfuctions. Inserting resolutions of the identity, J dq \q){q\ = 1 and making 
use of the closure relation for mapping states, X^a(9'I"^a) ("IaI?) = S{q — q'), the correlation function 
takes the form 




\X'titj.'=l 




(14) 



CABit) 




(15) 
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It is important to note that the dimensionahty of the mapping coordinate space representation is fixed 
by the number of quantum states. 

A phase space description of this coordinate representation can be obtained by introducing Wigner 
transforms of the matrix elements. Making the change of variables, qi = ri — zi/2, q2 = ri+ zi/2, etc., 
and, using the analog of Eq. ^ for the mapping coordinates, we find 

CAB{t) = (2TTnY^b+N) j dXidX2dxidX2 {peqA)wmiXl, Xi) 

xBwmiX2,X2)W{Xi,X2,Xi,X2,t), (16) 

where, in analogy with the bath phase space terminology, x = {r,p). The full Wigner representation of 
the mapping spectral density, W, is defined as 

W{Xi,X2,xi,X2,t) = j dZidZ2dzidz2 in + yKi?i + |i|e*^™*/^ji?2 - Y)|r2 - y) 

X {r2 + y |(i?2 + ^|e-*^™*/^|i?i - ^)|ri - |.)e-i(^-^i+^2-^2+Pi-i+f2-2). (17) 

All of these manipulations have served simply to cast the exact expression for the quantum correlation 
function into an equivalent exact form involving a phase-space-like representation for both the bath and 
the mapping version of the quantum subsystem degrees of freedom. There are advantages to defining W 
and deriving its dynamics rather than directly considering the dynamics of the correlation function. Not 
only are the algebraic manipulations simplified, but also all correlation functions for a specific system 
share the same spectral density. This fully quantum problem is still intractable for complex many-body 
systems. We now turn to the evaluation of the spectral density function in the quantum-classical limit 
that will provide the basis for simulation algorithms of the dynamics. 



5 Quantum- Classical Dynamics for W 

Given the development presented above, the dynamical problem consists in finding a tractable equation 
of motion for the spectral density function W{Xi,X2,xi,X2,t) and constructing an algorithm for its 
solution. The equation of motion can be derived by differentiating W with respect to time. However, 
to derive the quantum-classical equation for W we must introduce approximations. For this purpose we 
make use of an expansion in the small mass ratio = (m/M)^/^, where it is assumed that the charac- 
teristic mass M of the bath particles is much larger than m for the quantum subsystem particles. The 
manner in which this expansion is carried out is analogous to that discussed earlier for quantum Liou- 
ville equation [3j. More specifically, given the energy eo, time to = ^/^o, and length Am = {'h'^ /meo)^^"^ 
units, we scale the light and heavy particle momenta, respectively, with pm = (mXm/to) = (meo)^^'^ and 
Pm = {Meo)^^"^. Note that the only difference between subsystem and bath particles is scaling their 
momenta with different factors. As the subsystem and bath are in thermal equilibrium their average 
kinetic energies are equal and therefore, on average, p/P = /x, so that after scaling both subsystem and 
bath momenta have the same order of magnitude. After scaling variables, we obtain 

W'{R[,Pi,R'„Pi,r[,p[,r'„p',;t') = J dZ[dZ'^dz'M{r[ + + ^W^'-''\R'^ - ^)\r', - |) 

where a prime means that the variable is divided by the corresponding scaling factor. To avoid cum- 
bersome notation, we drop the primes in the following relations. Differentiation of the scaled form of 
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W yields the equation of motion, 



dW{t) 



i / dZidZ2dzidz2dQdq 



{ri + ^\{Ri+^\Hm\Q)\q) 



2 



-in + y|(i?i + y |e^^-*|i?2 - Y)|r2 - f )(r2 + y |(i?2 + f|e-^^-*|g)|g) 



(19) 



The mapping Hamiltonian (Eq. (jl3l) ) in scaled variables is 



+ Vb{R) + Exx' 



r\r\' 



PxPx' ~ <^AA')) while the scaled quantum mechanical momentum operators are P = ^-§q and p = —i-§^- 
Substituting the Hamiltonian and momentum operators into Eq. p9|) yields, 



dW{t) 
dt 



/ dZidZ2dQdzidz2dqe~'^^'^^+^^^^'>e-'^P'-''+P'-^^^ 



92 



1 d d 



xS{Q -Ri- ^)5{q - n - ^) \{q\{Q\e^"-'\R2 - ^)\r2 - ^) 



£2^ 
2 



x(r2 + yKi?2+ 2 



2 



AA' 



dq\i 



xS{Q - R, + ^)6{q - n + |)}(n + + ^|e^^-'|i22 - ^)\r2 - |) 



^2 , 



x(^2 + yKii2+ 2 



(20) 



where the change of variable Zi = Z\l [i was performed in order to transfer the /i dependence from the 
exponential to the argument of the potential energy terms, after integrations over Q and q have been 
performed. It is more convenient to expand a potential energy term around a small argument rather 
than to deal with an oscillatory exponential. 

We begin the calculation by performing the integrals over q and Q, followed by Taylor expansion 
of the potential energy terms around = 0, keeping terms up to the first order in /i. Finally, we 
introduce the definition of W into the resulting expression. The algebra is lengthy so it is presented 
in the Appendix. While the algebra leading to the result is lengthy, the final equation of motion is 
relatively simple. In unsealed coordinates it takes the form. 



dW{t) 



Y^hxx'iRi 



d 



d 



rix 



XX' 



8 ^ 

AA' 

= iCm{xi,Xi)W{t). 



dpix' 
d 



dRi 



drix 

d d 
+ 



d 



d 



d 



drix drix' dpiy dpix ) dP] 



M dRi 
d 

-Wit) 



dRi dPi 



)W{t) 



(21) 
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Here the Wigner transform of the Hamiltonian (jl3p is given by 

^"^ = + ^^^^^ + ^ ^^^fr^^'""'^^' " ^^^^ 

AA' 

This equation of motion is one of the principal results of this paper. It results from an expansion 
of the evolution operator for W to order fj, and is equivalent to quantum-classical Liouville dynamics 
for the spectral density function [T7]. The first term in Eq. (j2ip represents the subsystem dynamics 
of the spectral density in the mapping phase space. The second term is the dynamics of the spectral 
density due to classical evolution of the bath degrees of freedom and the third term is a higher order 
correlation between the quantum mapping and classical degrees of freedom. The last equality defines 
the quantum-classical Liouville operator in the mapping basis. This equation must be solved subject 
to the initial condition, 

WiO) = {2T:n)^^^+^h{ri-r2)6{Ri-R2)6{pi-p2)5{Pi-P2). (23) 

In Eq. (jl9p . since Hm commutes with propagator e*^™*, when differentiating with respect to time, 
Hm could be placed on either side of the propagator. In this derivation, we chose to put it to the left 
of the propagator in the first term in Eq. (jl9p and to the right of the propagator in the second term. If 
instead one places Hm to the right of the propagator in the first term and to the left of the propagator 
in the second term we can obtain an alternate form of the equation of motiorH. The manipulations are 
similar to those described above and in the Appendix and are not repeated her. The resulting equation 
of motion is 

= +^E^AA'(i?2)hA^-m^lw^(t)-(^-^-7^-^)H^(i) 



AA' 



^S^ dhxy{R2) /_5 d d d \ d 

8 ^ dR2 ■ \dr2x dr2X' dp2X' dp2x)dP2 ^ ^ 

= -iCm{x2,X2)W{t). (24) 

This alternate form of the dynamics not only shows the symmetry of the dynamics of W in its variables 
but also allows us to move the time evolution from W to the observable in the correlation function 
expression. The formal solution of Eq. (I24p is 

W{Xi,X2,xi,X2,t) = e-'^™("2'^2)*iy(Xi,X2,xi,X2,0). (25) 

These results will be used in the next section to derive a quantum-classical approximation to the time 
correlation function. 



6 Quantum- Classical Correlation Function 

We can now employ the results of the last two sections to find an expression for the quantum-classical 
approximation to the correlation function in the mapping basis. Using Eq. (125p in Eq. ()16p . we have 

CABit) = ^27rhy^^b+N) j dXidX2dxidX2 {peqA)wmiXl,Xi) 

xBwm{X2, a;2)e-^^-("2'^2)*I^(Xi, X2, XI, X2, 0). (26) 

■^If firn is placed in other locations, say, either to the right or left in both terms, the evolution operator for W is the 
mean of the two forms discussed in the text. 
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Performing an integration by parts to move the evolution operator onto the and integrating over 

the coordinates with subscript 1, making use of the delta functions in the initial value of W, yields, 

CAsit) = j dXdx {peqA)wmiX,x)BwmiX,X,t), (27) 

where Bwm{^, x, t) = e'-^"'^^'^^^ Bwm{X , x). (We have dropped the subscripts on the phase space coor- 
dinates since this notation is no longer needed.) This expression contains the full quantum equilibrium 
structure of the subsystem and bath and the quantum-classical Liouville evolution of the operator B 
in the mapping representation. 

The evolution equation for By[rm{X,x,t) that one obtains in this derivation is identical to that 
found earlier in the calculation of the average value of an observable |10j . This evolution equation can 
be written in the form 

-J^Bwrnix, X,t) = iCrnBwm{t) = -{Hm, Bwmit)}x,X (28) 

h ^ dhxx' ^ ^ , ^ ^^Ar (f\ 

+ 8^ dR '^drx'drx^ dpx'dpx^dP "^"^^ >' 

AA 

where { , }x,x denotes a Poisson bracket in the full mapping-bath phase space of the system. The 
quantum-classical Liouville operator in the mapping basis can be decomposed into two terms, iCm = 
iC^ + iC'^ where 



iCl, = -{Hm, }x,X, (29) 
^ m - fil^ Qji '^drx'drx dpx'dpx'dP' 



XX' 

The iC^ evolution operator, corresponding to the Poisson bracket in Eq. (p8]) . leads to a classical- like 
evolution of the coupled dynamics of the quantum mapping and classical bath phase space variables 
that can be simulated by Newtonian trajectories. The force field that the classical variables feel varies 
with time as a result of the dependence of the forces on the mapping phase space variables. If iC'^ 
is dropped in the evolution equation, Bwrnit) has a solution in terms of characteristics. The set of 
ordinary differential equations that determines its solution is 



drxit) 



dt h 

X 



i J^/lAA'(i?(0)PA'(t), 



dpx{t) 



ij;/iAA'(i?(t))rA'(i), 



dt n 

A' 

dm ^ P(t) dm^_dHrn_ . . 

dt M ' dt dR{ty ^ ' 

The utility of this approximation to Eq. (I28p is a topic of current research. Tests of its accuracy have 
been carried out on the spin-boson model where Eq. (j28p is equivalent to full quantum dynamics and are 
being carried out on other model systems with nonlinear coupling between the quantum and classical 
degrees of freedom where Eq. ()28p is not exact. For the spin-boson model Eq. (j30p yields results that 
are indistinguishable from the known exact quantum results for this system |10j . Tests being carried 
out on other model systems show that while the results are often in close accord with exact quantum 
results, sometimes discrepancies are observed that cannot be ascribed to approximations in Eq. (I28p 
for these more general interactions and must be attributed to the use of Eq. ()30p for the dynamics. 
Consequently, further research is underway to fully characterize the contributions arising from iC'^ and 
construct algorithms that account for its action. 
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7 Conclusion 



In the adiabatic basis the correlation function may be evaluated in terms of an ensemble of surface- 
hopping trajectories O |5]. While this simulation scheme has been used in a number of applications 
it is difficult to obtain accurate results for long times due the presence of oscillating terms in the 
Monte Carlo sampling. In this paper, we show that it is possible to reformulate the calculation of the 
correlation function in the mapping basis. Simulation schemes based on this formulation do not suffer 
from some of the problems that arise in the implementation of surface- hopping schemes. 

In order to evaluate the expression for the correlation function obtained in this paper two ingredients 
are required. The expression involves an average over the full quantum equilibrium structure. As briefly 
discussed above, the computation of quantum equilibrium structure is a more tractable problem than 
the calculation of quantum dynamics. Nevertheless, approximations are usually required to evaluate 
the equilibrium structure. In the high temperature limit it is possible to derive analytical expressions 
that are useful in many applications [321 [33| [34] . 

The main result of this paper is the quantum-classical expression for the correlation function that 
involves time evolution of the quantum subsystem in the mapping basis. Provided the quantum- 
classical evolution is approximated by iCm — i^m the time evolution of the dynamical variable in 
the correlation function can be computed easily by solving a set of Newtonian-like equations. Thus, 
difficulties associated with the accumulation of Monte Carlo weights in the evaluation of an oscillatory 
function that arise in the surface-hopping solution of the quantum-classical Liouville equation in the 
adiabatic basis are by-passed. Of course, this simple scheme relies on the ability to neglect iC^, which 
accounts for higher order correlations in the dynamics. Thus, the focus of future research is on the 
characterization of the nature of the dynamics generated by i>C^ and the construction of simulation 
algorithms that account for its presence. The results in this paper form the basis for future applications 
to the calculation of transport properties, such as rate constants for nonadiabtic chemical reactions. 
Acknowledgement: This work was in part supported by a grant from the Natural Science and Engineer- 
ing Research Council of Canada. AN acknowledges the support from the Lachlan Gilchrist fellowship. 



Appendix: Derivation of the Dynamics 

In this Appendix, we give the details of the calculations needed to obtain the quantum-classical evolution 
equation for the spectral density function. There are ten terms in Eq. (I20p and, because of the symmetry 
of the expression, it is convenient to group term i with term i-|-5 and evaluate the contributions group 
by group. First, we perform the integrals over Q and q. The 1st and 6th terms contain and 

respectively, after this integration. Integration by parts with respect to Zi and Z2 and summation 
of the results yield /iPi • The 2nd and 7th terms involve bath potentials. After expansion of 

Vb{Ri ± in the first and third terms in the series, which are proportional to and cancel 
and the second terms, which are proportional to /i, yield " f^- The contribution from the 

5th and 10th terms derived similarly to yield ^ ^^r^ ' f^- The derivations of the third and eighth 
and, also, forth and fifth groups of terms are more complicated and are presented in the following two 
subsections. 
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3rd and 8th terms: 

The sum of the 3rd and 8th contributions to the time derivative of W is 



/dW\ 

Km) 



3+8 



hxx'{Ri) + 
h\\'{Ri 



/ dZidZ2dzidz2e-'^^^-^^+^^-^^+P^-^^+P^- 
XX' •' 

dhxx'{Ri) I^Zi 



22j 



dRi 2 
dRi ' 2 - 



nxrix' + rixzix'/2 + zixrix'/2 + zu^u'/^ 



rixrix' - rixzix'/2 - zixrix'/2 + zu^u'/^ 



x(ri + |Ki^i + ^|e^^™*|i?.-^>h-|> 



x(r2 + ||(i?. + ^|e-Mi^,-^)h-|), 

where we have performed a McLauren expansion of hxx' in the small parameter f^, and retained the 
first two terms. Finally, using the fact that a partial differential with respect to momentum acting on 
the exponential term has the same effect as multiplication by the variable zi or Zi that results from 
the expansion, we have 



\ Ot /3+8 



IfJ, 



2Nt 



I dZidZ2dzidz2^ihxx'{Ri)[rix-Q^ + rix 



/X dhxx' {Ri 



d 
'P: 

5 1 d d 



d 



dp 



IX 



2 dRi dPi L 2dpixdpix 



-i{Pl-Zi+P2-Z2+pi-Zl+P2-Z2) 



x{ri+"-^\{Ri + ^\e'"-' 



\R2-—)\r2--) 



AA' 



dp 
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lA' 



dp 



\W 



lA' 



fj. sr^ dhxx'{Ri) r„ 1 9 9 n 5 



AA' 



2 9piA 5piA' 9Pl 



where in the last equality the expression for the scaled W is inserted. 
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4th and 9th terms: 

The sum of the 4th and 9th terms in Eq. (j20p is 



/dW\ 

\~dr) 



2Nt 



4+9 



dZidZ2dzidz2 < |^/iaa'(-^ 

AA' I 



1 + 



/iZi d d 



x(n + ||(i?i + ^KMi?,-^)|r,-|) 
x(r. + ||(i?. + ^|e---*|i?,-4^)|n-|; 



d 



2 ^9(ri - f )a 9(ri - f ); 



x(ri + ||(i?i + ^|e^^-*| 



2 /A i 2 



> ki 

2 ' 2 ' 



_R2 - — )|r2 - — ) L-*(^l-^l+'P2-^2+Pl-21+P2-22) 



2Nt 



/ (iZi(iZ2(izi(i22< (/iAA'(-^i) + 

AA' I 



dRi 2 



) 



x(r2 + ||(i?2 + ^|e---'|i?i-4^)|n-|) 



{hxx'{Ri) 
d d 



dhxx'{Ri)i^Z 



dRi 

^2, 



3 ^ ■(.2 + f|{fl. + i^|.-"»"'lft-^)ln , 

orix orix' 2 2 2 2 



fiZi 



(31) 



Again in the second equahty we have carried out a McLauren expansion of hxx' to first order in /i. 
Equation (j3ip itself has four terms. In the sum of the 1st and 3rd subcontributions one of the partial 
derivatives over r may be replaced by a partial derivative over z and an integration by parts may be 
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carried out. We find 



2Nb 



AA' 



x(^-.m)(r2 + ||(^2 + — le 



1-— )h-y, 



5 



(32) 



AA' 



Here we see that the two terms with a momentum multipher are added together while the other terms 
cancel each other. 

The 2nd and 4th sub contributions in Eq. (|31|) are treated in the following way: We replace multi- 
plication by the variable Z with a partial differentiation with respect to the momentum. Furthermore, 
the sum of these contributions is written as one half the sum of two equal contributions, the expressions 
appearing in Eq. (jSip and the same expression but with the differentials over r replaced with those 
over z. Thus, the sum of the 2nd and 4th contributions can be written as 



AA' 

d d 



dRi 



dZidZ2dzidz2ii^)e-'^^'-^'+^'-^'+P'''+P^-'''^ 
oPi 



^2 , 



^lZ2. 



x(r2 + y|(i?2 + — le 



\R 



Zl. 



+ (n + ||(i?i + 4^|e--'|i2,-^)|r,-|) 



+4 



_d_ d_ 

drix drix 
d d 



R^ 



^Z 



1 V I Zl 



^^(n + tl(i?i + ^!e^^-*|i^2-^ 
ozi\ ozix' 2 2 2 



£l 
2 

d d 



I1Z2 , I ^2 V 

2 )N--) 



\l ^2v 

)k2-y) 



(33) 



If in both expressions we replace ^§^B + with ^^^^ - 1^1^ - 1^ If > the cross terms cancel. 



dxdy 



dxdy dx dy dy dx 



13 



Finally, integrating by parts over zix and zw , we get 

dRi '^dnxdnx' ^'^^-'^'^'dPi 



jj. sr^ dhx\i{Ri) . d d ^ N 5 TT^ /o.^ 



AA' 

Summing all the contributions from the above calculations, we find the comparatively simple result 



dW 
~~dt 



d d 1 

opix' orix' J 



-Y,hxx'{Ri) 

AA' 

_^ r, dW dVB{Rl) dW f,^ dhxx'{Rl) , ^ A ^ ^ M/ 

u^dhxy{Ri)_ d_ d d_\_d_ 

^ 8 ^ dR, '\dnx dnx' dpix dpix' ) dP^ ' ^ ' 

Restoring unsealed coordinates we have Eq. (pT]) in the text. 
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